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Abstract 

Sum-rate maximization in two-way amplify-and-forward (AF) multiple-input multiple-output (MIMO) 
relaying belongs to the class of difference-of-convex functions (DC) programming problems. DC pro- 
gramming problems occur as well in other signal processing applications and are typically solved 
using different modifications of the branch-and-bound method. This method, however, does not have 
any polynomial time complexity guarantees. In this paper, we show that a class of DC programming 
problems, to which the sum-rate maximization in two-way MIMO relaying belongs, can be solved very 
efficiently in polynomial time, and develop two algorithms. The objective function of the problem is 
represented as a product of quadratic ratios and parameterized so that its convex part (versus the concave 
part) contains only one (or two) optimization variables. One of the algorithms is called POlynomial- 
Time DC (POTDC) and is based on semi-definite programming (SDP) relaxation, linearization, and an 
iterative search over a single parameter. The other algorithm is called RAte-maximization via Generalized 
Eigenvectors (RAGES) and is based on the generalized eigenvectors method and an iterative search over 
two (or one, in its approximate version) optimization variables. We also derive an upper-bound for the 
optimal values of the corresponding optimization problem and show by simulations that this upper-bound 
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can be achieved by both algorithms. The proposed methods for maximizing the sum-rate in the two-way 
AF MIMO relaying system are shown to be superior to other state-of-the-art algorithms. 

Index Terms 

Difference of convex functions programming, Non-convex programming, Semi-definite programming 
relaxation, Sum-rate maximization, Two way relaying 

I. Introduction 

Two-way relaying has recently attracted a significant research interest due to its ability to overcome the 
drawback of conventional one-way relaying, that is, the factor of 1/2 loss in the rate (H, |2). Moreover, 
two-way relaying can be viewed as a certain form of network coding Q which allows to reduce the 
number of time slots used for the transmission in one-way relaying by relaxing the requirement of 
'orthogonal/non-interfering' transmissions between the terminals and the relay |4!]. Specifically, simulta- 
neous transmissions by the terminals to the relay on the same frequencies are allowed in the first time 
slot, while a combined signal is broadcasted by the relay in the second time slot. In contract to the one- 
way relaying case, the rate-optimal strategy for two-way relaying is in general unknown 0. However, 
some efficient strategies have been developed. Depending on the ability of the relay to regenerate/decode 
the signals from the terminals, several two-way transmission protocols have been introduced and studied. 
The regenerative relay adopts the decode-and-forward protocol and performs the decoding process at the 
relay [6], while the non-regenerative relay typically adopts a form of amplify-and-forward (AF) protocol 
and does not perform decoding at the relay, but amplifies and possibly beamforms or precodes the signals 
to retransmit them back to the terminals 10, Q, JH. The advantages of the latter are a smaller delay in 
the transmission and lower hardware complexity of the relay. 

In this paper, we consider the AF two-way relaying system with two terminals equipped with a single 
antenna and one relay with multiple antennas. The task is to find the relay transmit strategy that maximizes 
the sum rate of both terminals. This is a basic model which can be extended in many ways. The significant 
advantage of considering this basic model is that the corresponding capacity region is discussed in the 
existing literature in O. It enables us to concentrate on the mathematical issues of the corresponding 
optimization problem which are of significant and ubiquitous interest. 

We show that the optimization problem of finding the relay amplification matrix for the considered AF 
two-way relaying system is equivalent to finding the maximum of the product of quadratic ratios under 
a quadratic power constraint on the available power at the relay. Such a problem belongs to the class 



February 14, 2012 



DRAFT 



3 



of the so-called difference-of-convex functions (DC) programming problems. It is worth stressing that 
DC programming problems are very common in signal processing and, in particular, signal processing 
for communications. For example, the robust adaptive beamforming for the general-rank (distributed 
source) signal model with a positive semi-definite constraint can be shown to belong to the class of DC 
programming problems O, |[T0l . Specifically, the constraint in the corresponding optimization problem 
is the difference of two weighted norm functions. The power control for wireless cellular systems is 
also a DC programming problem when the the rate is used as a utility function ifTTTl . Similarly, the 
dynamic spectrum management for digital subscriber lines lfT2l as well as the problems of finding 
the weighted sum-rate point, the proportional-fairness operating point, and the max-min optimal point 
(egalitarian solution) for the two-user multiple input single output (MISO) interference channel |[T3l are 
all DC programming problems. The typical approach for solving such problems is the use of various 
modifications of the branch-and-bound method |[T3l - |[T9l that is an efficient global optimization method. 
The branch-and-bound method is known to work well especially for the case of monotonic functions, 
i.e., the case which is typically encountered in signal processing and, in particular, signal processing for 
communications. However, it does not have any worst-case polynomial complexity guarantees, which 
significantly limits or essentially prohibits its applicability in practical communication systems. Thus, 
methods with guaranteed polynomial-time complexity that can solve different types of DC programming 
problems are of a fundamental importance. 

In the last decade, a significant progress has occurred in the application of optimization theory in 
signal processing and communications. Some of those results are relevant for the considered problem 
of maximizing constrained product of quadratic ratios E01 - E31 . The worst-case-based robust adaptive 
beamforming problem is known to belong to the class of second-order cone (SOC) programming problems 
EU1 largely due to the fact that the output signal-plus-interference-to-noise ratio (SINR) of adaptive 
beamforming is unchanged when the beamforming vector undergoes an arbitrary phase rotation. This 
allows to simplify the single worst-case distortionless response constraint of the optimization problem 
into the form of a SOC constraint. The situation is significantly more complicated in the case of multiple 
constraints of the same type as the constraint in EU1 when a single rotation of the beamforming vector 
is not sufficient to satisfy all constraints simultaneously. This situation is successfully addressed in Ell 
by considering the semi-definite programming (SDP) relaxation technique. The SDP relaxation technique 
has been then further developed and studied in, for example, E2l . E3l and other works. Interestingly, 
the work E3l considers the fractional quadratically constrained quadratic programming (QCQP) problem 
that is closest to the one addressed in this paper with the significant difference though that the objective 
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in |[23l contains only a single quadratic ratio that simplifies the problem dramatically. 

In this paper, we develop polynomial time algorithms for finding the globally optimal solution of 
a class of non-convex DC programming problems, e.g., the maximization of a product of quadratic 
ratios under a quadratic constraint. This problem precisely corresponds to the sum-rate maximization in 
two-way AF MIMO relaying. Our algorithms use such parameterizations of the objective function that 
its convex part (versus the concave part) contains only one (or two) optimization variables. One of the 
proposed algorithms is named POlynomial-Time DC (POTDC) and is based on semi-definite programming 
(SDP) relaxation, linearization, and an iterative search over a single parameter^ The POTDC algorithm 
is rigorous and finds the global maximum of the considered problem. Indeed, the solution given by this 
algorithm coincides with the newly developed upper-bound for the optimal value of the problem. The 
other algorithm is called RAte-maximization via Generalized Eigenvectors (RAGES) and is based on 
the generalized eigenvectors method and an iterative search over two (or one, in its approximate version) 
optimization variables^ The RAGES algorithm is somewhat heuristic in its approximate version, but may 
enjoy a lower complexity. 

The rest of the paper is organized as follows. The two-way AF MIMO relaying system model is given 
in Section II while the sum-rate optimization problem for the corresponding system is formulated in 
Section III. The POTDC algorithm for solving the corresponding sum-rate maximization is developed in 
Section IV and an upper-bound for the optimal value of the maximization problem is found in Section V. 
In Section VI, the RAGES algorithm is developed and investigated. Simulation results are reported in 
Section VII followed by the conclusions. This paper is reproducible research ||26l and the software needed 
to generate the simulation results will be provided to the IEEE Xplore together with the paper upon its 
acceptance. 

II. System Model 

We consider a two-way relaying system with two single-antenna terminals and an amplify-and-forward 
(AF) relay equipped with Mr antennas. Fig. 1 shows the system we study in the paper. In the first 
transmission phase, both terminals transmit to the relay. Assuming frequency-flat quasi-static block fading, 
the received signal at the relay can be expressed as 

r = hj /) • xi + h { 2 f) -X2 + I1R (1) 

'Some preliminary results on the POTDC algorithm have been submitted to ICASSP'12 1241 . 
2 Some preliminary results on the RAGES algorithm have been presented in 1251 . 
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where = [/i^i, . . . , hi t M R ] T £ C RIr represents the (forward) channel vector between terminal i and 
the relay, Xi is the transmitted symbol from terminal i, iir € C Mr denotes the additive noise component 
at the relay, and (-) T stands for the transpose of a vector or a matrix. Let Pp^ = E{|xj| 2 } be the average 
transmit power of terminal i, Hn,r = E{riR ■ n|[} be the noise covariance matrix at the relay, E{-} 
denoting the mathematical expectation, and (-) H standing for the Hermitian transpose of a vector or a 
matrix. For the special case of white noise we have Hn,r = Pn,r • Ia/ r where Pn,r = tr(R/v,#)/MR 
and Im r is the identity matrix of size Mr x Mr. 

The relay amplifies the received signal by multiplying it with a relay amplification matrix G € 
C sX s , i.e., it transmits the signal r = G • r. The transmit power used by the relay can be expressed 
as 

E{||f||l} = E{tr{G-r-r^-G H }} 

= tr{G-R /? -G // } =ti{G H -G-YLr) (2) 

where || • ||2 denotes the Euclidian norm of a vector and Rr = E{r • r^} is the covariance matrix of r 
which is given by 

R R = h[ f) ■ (h^) H ■ P T)1 + h ( 2 f) ■ (h ( 2 f) Y ■ P T , 2 + R N)R . (3) 

Next, we use the equality tr(A H ■ B) = vec(A) H • vec(B), which holds for any arbitrary square matrices 
A and B, and where vec(-) stands for the vectorization operation that transforms a matrix into a long 
vector stacking the columns of the matrix one after another. Then, the total transmit power of the relay 
© can be equivalently expressed as 

^{llrlli} = v ec(G) ff -vec(G-R R ). (4) 

Finally, using the equality vec(A • B) = (B T ® I/,) • vec(A), which is valid for any arbitrary square 
matrices Al x l and T$lxL, and where <8> denotes the Kronecker product, © can be equivalently rewritten 
as the following quadratic form 

E{||r|||} = g H ■ (R£ I Mr ) g = • Q • g (5) 
Q 

where g = vec{G}. 

In the second phase, the terminals receive the relay's transmission via the (backward) channels (h^) T 
and (h^) T (in the special case when reciprocity holds we have b!p = hf' for i = 1, 2). Consequently, 
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the received signals yi, i = 1, 2 at both terminals can be expressed, respectively, as 

2/2 



h^-x 1 + h^-x 2 + h 1 

^2 2 ' X 2 + «2 1 • x l + n 2 



(6) 
(V) 



where h 



and hi 



(e) 



hf») T 



T 



( f) 

■ G ■ h - is the effective channel between terminal i and terminal j for i, j = 1, 2 



• G • n^j + nj represents the effective noise contribution at terminal i which comprises 
the terminal's own noise as well as the forwarded relay noise. The first term in the received signal of 
each terminal represents the self-interference, which can be subtracted by the terminal since its own 
transmitted signal is known. The required channel knowledge for this step can be easily obtained, for 
example, via the Least Squares (LS) compound channel estimator described in ||27l . 

After the cancellation of the self-interference, the two-way relaying system is decoupled into two 
parallel single-user SISO systems. Consequently, the rate T{ of terminal i can be expressed as 



1 



Id 1 + 



1 



Id 



P 



(8) 



N., 



where ld(-) denotes the logarithm of base two, Pr^ and Pjv,j are the powers of the desired signal 
and the effective noise term at terminal i, respectively, and Pm = Pm + Pn i- Specifically, Pjh = 



E 



, Pr- 



E 



h!£\ ■ xi k and P N)i = E { |nj| 2 } for i = 1,2. Note that the factor 1/2 



results from the two time slots needed for the bidirectional transmission. The powers of the desired signal 
and the effective noise term at terminal i can be equivalently expressed as 




G h 



(/) 



• G • h\ 



(/) 



(9) 
(10) 

(11) 
(12) 



where the expectation is taken with respect to the transmit signals and also the additional noise terms. 
Moreover, these powers can be further expressed as quadratic forms in g. For this goal, first note that 
by using the following equality 

vec(A • B • C) = (C T <g> A) • vec(B) (13) 
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which is valid for any arbitrary matrices A, B and C of compatible dimensions, the term fhj- 6 ^ G h 
can be modified as follows 



(/) 

3 



h< 6) Y 



Gh 



(/) 



j 



vec ( h 



• G • 



h 



(/) 



vec(G). 



Using ( fT31 ), the power of the desired signal at the first terminal can be expressed as 



PR, 



.(/) 



h?>) 



T 



H 



h 



(/) 



.(*>) 



Pt;. 



(14) 
(15) 

(16) 



Finally, applying the equahty (A <g> B) • (C <g) D) = (A • C) <g> (B • D) to (fl6l ) which is vahd for any 
arbitrary matrices A, B, C and D of agreed dimensions, Pr : \ can be expressed as the following quadratic 
form 



Pr,i = E H ■ 



h 



(/) 



h 



(/) 



H 



I h? ■ (h? 



11 



T 



g • Pr- 



ill) 



Similarly, Pr^ can be obtained. 

By defining the matrices K^,!, K12 as follows 



K 



2.1 



h 



(/) 



K 



1,2 



h 



(/) 



h 



(/) 



H 



H 



h<" 



h 



(6) 



h 



?>)" 



the powers of the desired signal can be expressed as 

Pr,i = • K 2 ,i • g • P T ,2 
P R , 2 = g H ■Ki.a-g-Pr.i. 



(18) 
(19) 

(20) 
(21) 



As the last step, the effective noise Pjvi can be converted into a quadratic form through the following 
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train of equalities 



N,i 



E 



h 



(6) 



T 



G ■ n R + n, 



G • K NtR ■ G H ( hf r ) 



tr 



R 



iV.if, 



Pi 



vec(G)^ • vec (hj 



.(6) 



(6) 



' G • ~Rn,R ] + PN,i 
T 



vec(G)" • (r n>r ® (hf) (h| 6) )*) • vec(G) + Pat, 



(22) 
(23) 

(24) 
(25) 



where d23l is obtained from (1221) by applying the the equality tr(A H • B) = vec(A)^ • vec(B), which 
is valid for any arbitrary square matrices matrices AlxL and HlxL, equation (1241 is obtained from (l23l 
by applying the equality (fl"3l . and the matrix Jj is defined as 



R 



AT/? 



(26) 



III. Problem statement 



Our goal is to find the relay amplification matrix G which maximizes the sum-rate r\ + r 2 subject 
to a power constraint at the relay. For convenience we express the objective function and its solution in 
terms of g = vec{G}. Then the power constrained sum-rate maximization problem can be expressed as 



gopt = arg max ( n + r 2 ) subject to g H ■ Q • g < P t ,r 



(27) 



where Pt,r is the allowed transmit power at the relay. Using the definitions from the previous section, 
this optimization problem can be rewritten as 



1 



gopt = arg max -Id 

g|g H -Qg<^T,fi 1 



arg max [ 1 + 
g|g H -Q g<^T,H \ Pn,i 



Pn,1 J \ P/V,2 / 



Pi 



N,2 



arg max 



Pra Pr,2 



(28) 



(29) 



s\s h -Q s<Pt,r Pv,1 Pv,2 

where we have used the fact that 0.5 • ld(x) is a monotonic function in x € M + and Pr^% = 1,2 is 
defined after ([8]>. 
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It is worth noting that the inequality constraint in this optimization problem has to be active at the 
optimal point. This can be easily shown by contradiction. Assume g op t satisfies g^p t • Q • g op t < Pt,r- 
Then we can find a constant c > 1 such that g opt = c • g opt satisfies g^ t • Q • g opt = Pt,r- However, 
inserting g opt in the objective function of d28l) . we obtain 



1 c 2 ■ gg, t K 2 ,ig opt Pr i2 \ f c 2 • g^pt K i,2go P t-Pr,i 



V c ~ • gSt J igo P t + Pn,i) \ c 2 • g# t J 2 g op t + Pn,2 ) 

( gg,t K 2,lgopt^T,2 \ ( gg,t K 1,2 gopt Pr,l \ 

+ _H t- . Pn,i j ' \ + „H T _ ,Ps.») { ' 

V g op t J Igopt + J \ g op t J 2gopt + / 

which is monotonically increasing in c. Since we have c > 1, the vector g opt provides a larger value of 
the objective functions than g opt which contradicts the assumption that g opt was optimal. 

As a result, we have shown that the optimal vector g opt must satisfy the total power constraint of 
the problem with equality, i.e., g^L • Q • g opt = Pt,r- Using this fact, the inequality constraint in the 
problem d29l ) can be replaced by the constraint g^ • Q • g = Pt,r- This enables us to substitute the 
constant term Pjy,i, which appears in the effective noise power at terminal i (l25l) . with the quadratic 
term of g^L • Q • g opt • (Pn,i/ 'Pt,r)- This leads to an equivalent homogeneous expression for the ratio 
of PR : i/PN,i,i = 1)2. Thus, by using such substitution, Pn,%, i = 1,2 from d25l) can be equivalently 
written as 

P Nyi = g H ■ B, • g, i = l,2 (31) 



where Bj is given by 



Bi = Ji + ^1 Q. (32) 

Pt,r 



Inserting (f20]), (1211 ). and 02l ) into d29l ), the optimization problem becomes 

g H ■ Ai • g g H ■ A 2 • g 

g op t = arg max y^—- (33) 

where we have defined the new matrices Ai = K 21 • Pti + Bi and A 2 = K 12 • Pt \ + B 2 . 

As a final simplifying step we observe that the objective function of (l33l) is homogeneous in g, meaning 
that an arbitrary rescaling of g has no effect on the value of the objective functions. Consequently, the 
equality constraint can be dropped completely as any solution to the unconstrained problem can be 
rescaled to meet the equality constraint without any loss in terms of the objective functions. Therefore, 
the final form of our problem statement is given by 

g H ■ Ai • g g^ • A 2 • 

Lax 

g 



go P t = arg max „ 3-— . (34) 

% H ■ Bi • g g H ■ B 2 • g 
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Note that from their definitions it is obvious that Aj, i = 1,2 and B^, i = 1,2 are positive definite 
matrices. Therefore, the optimization problem (l34l ) can be interpreted as the product of two Rayleigh 
quotients. Moreover, it can be expressed as a DC programming problem. Indeed, as we will show later 
in details, by expressing the problem (l34l) as a rank constrained problem and then dropping the rank 
constraint and also taking the logarithm of the objective function, the objective function of the resulting 
problem can be written as the summation of two concave functions with positive signs and two concave 
functions with negative signs. Thus, the objective of the equivalent problem is, in fact, the difference 
of convex functions which is in general non-convex, and the available algorithms in the literature for 
solving such DC programming problems are based on the so-called branch-and-bound method that does 
not have any polynomial time computational complexity guarantees lfT3l - |[T9ll . However, as we show next, 
the problem ( f34T > can be parameterized in such a way that there exist simple polynomial time solutions. 

IV. Polynomial-Time Solution for the Sum-Rate Maximization Problem in Two-Way 

AF MIMO Relaying 

Since the problem d34l) is homogenous, without loss of generality, we can fix the quadratic term 
g^ • Bi • g to be equal to one at the optimal point. By doing so and also by defining the additional 
variables r and (3, the problem (l34l) can be equivalently recast as 

max g H ■ Ai • g • 1 - 

e,T,P p 

g H • B x • g = 1 
g H • A 2 • g = r 

g H • B 2 • g = (35) 

Using the fact that the quadratic function g^ • Bi • g is set to one, one can easily check that the prob- 
lem 051 ) is feasible if and only if r € [A m j n (B 1 - 1 A 2 ),A inax (B^ 1 A2)] and (3 E [A^B^Bs), A max (Br ^2 
where A m ; n (-) and A max (-) denote the smallest and the largest eigenvalues operator, respectively. By 
introducing the matrix X = g-g^ and observing that for any arbitrary matrix Y, the equation g^ • Y-g = 



February 14, 2012 



DRAFT 



1 1 



tr(Y ■ g ■ g^) holds, the optimization problem (l35l ) can be equivalently expressed as 

7" 

max tr(Ai ■ X) • — 

X,r,/3 p 

tr(Bi • X) = 1 
tr(A 2 • X) = r 
tr(B 2 • X) = 

rank(X) = 1, X h 0. (36) 

In what follows, we explain the possibility of dropping the rank-one constraint in the problem (l36l) 
and then extracting the exact solution for the original problem (l36l ) based on the solution of the rank 
relaxed problem. To this end, let X T / 3 denote the optimal solution of the optimization problem (l36l ) with 
respect to X for fixed values of r and ft and without considering the rank-one constraint. It is known 
that the strong duality for a QCQP problem with three or less constraints is satisfied 11291 . Based on 
this fact, the strong duality holds for the problem (135T ). which for fixed variables r and ft is equivalent 
to QCQP with three constraints. Since the problem (l36l ) is equivalent to the problem (l35l ). the strong 
duality also holds for (1361 ) for fixed r and j3. As a result, a rank-one solution of the problem (1361 can 
always be constructed based on for fixed r and p. Thus, for fixed r and P, the optimal value 

of the problem (l36l ) with respect to X is independent of the rank-one constant. It enables us to drop 
the rank-one constraint in the problem (1361 ). solve the relaxed problem, and then construct an optimal 
rank-one solution once the optimal X opt , r opt , and p opt are obtained. Dropping the rank-one constraint 
results in the following optimization problem 

T 

max tr(Ai • X) • — 

X,r,/3 p 

tr(Bi • X) = 1 
tr(A 2 • X) = r 
tr(B 2 • X) = P 

X t 0. (37) 

Due to the fact that the matrix Ai is positive definite and X is positive semi-definite, the function 
tr(Ai • X) is always positive. The latter happens since the matrix X cannot be equal to a zero matrix 
due to the constraint tr(Bi • X) = 1. Moreover, since the values A m ; n (B| 1 ~ 1 A 2 ) and A m i n (B^ 1 B 2 ) are 
necessarily positive, the variables r and ft are also positive. The task of maximizing the objective function 
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in the problem (l37l) is equivalent to maximizing the logarithm of this objective function because log(rr) 
is a strictly increasing function and the objective function in (l37l) is positive. Therefore, the optimization 
problem 071 ) can be equivalently rewritten as 

max log(tr(Ai • X)) + log(r) — log(/3) 

X,t,/3 

tr(Bi • X) = 1 
tr(A 2 • X) = r 
tr(B 2 • X) = p 

X y (38) 

Note that dropping the rank-one constraint enabled us to write our optimization problem as a DC 
programming problem, where the fact that log(tr(Ai • X)) in the objective of (l38l) is a concave function 
is also considered. Although the problem d38l ) boils down to the known family of DC programming 
problems, still there exists no solution for such DC programming problems with guaranteed polynomial 
time complexity. However, the problem (l38l) has a very particular structure, such as, all the constraints 
are convex and the terms log(tr(Ai • X)) and log(r) in the objective are concave. Thus, the only term 
that makes the problem overall non-convex is the term — log(/3) in the objective. If — log(/3) is piece- 
wise linearized over a finite number of interval^, then the objective function becomes concave on these 
intervals and the whole problem (I38T ) becomes convex. The resulting convex problems over different 
linearization intervals for — log(/3) can be solved efficiently in polynomial time, and then, the suboptimal 
solution of the problem (l38l) can be found. The fact that such a solution is suboptimal follows from the 
linearization, which has a finite accuracy. The smaller the intervals are, the more accurate becomes the 
solution of (l38l) . This solution is also not the most efficient in terms of complexity. Thus, we develop 
another method (the POTDC algorithm) which makes it possible to solve the problem (I38T ) in a more 
efficient way. 

To fulfil this goal, we introduce a new additional variable t, which makes it possible to express the 

3 As explained before, the parameter /3 can take values only in a finite interval. Thus, a finite number of linearization intervals 
for — log(/3) is needed. 
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problem (I38T ) equivalently as 

max log(tr(Ai ■ X)) + log(r) - t 

X,r,/3 



tr(Bi • X) 


= 1 


tr(A 2 • X) 


= r 


tr(B 2 • X) 


= P 


log(/3) < t 




X t o. 





(39) 



The objective function of the optimization problem d39l) is concave and all the constraints except the 
constraint log(/3) < t are convex. Thus, we can develop an iterative method that is different to the 
aforementioned piece-wise linearization-based method, and is based on linearizing the non-convex term 
log(/3) in the constraint log(/3) < t around a suitably selected point in each iteration. More specifically, 
the linearizing point in each iteration is selected such that the iterative algorithm gets closer to optimal 
point in every iteration. Roughly speaking, the main idea of this iterative method is similar to the 
gradient based methods. In the first iteration, we start with an arbitrary point selected in the interval 
[A m j n (B^ 1 B2), A max (Bj" 1 B2)] and denoted as j3 c . Then the non-convex function log(/3) can be replaced 
by its linear approximation around this point /3 C , that is, 

log(/3) « logG8 c ) + ^-(P- Pc) (40) 

Pc 

which results in the following convex optimization problem 

max log(tr(Ai ■ X)) + log(r) - t 



tr(Bx 


•X) 


= 1 


tr(A 2 


•X) 


= r 


tr(B 2 


•X) 


= /3 



log(/3 c ) + i(/3-/? c )<t 

Pc 

X h 0. (41) 

The problem (14TT ) can be efficiently solved by means of the interior-point based numerical methods. 
Once the optimal solution of this problem in the first iteration, denoted as X^, and /3opt> i s f° un d, 
the algorithm proceeds to the second iteration by replacing the function log(/3) by its linear approximation 
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around /3^ t found from the previous (first) iteration. Fig. [2] shows how log(/3) is replaced by its linear 
approximation around j3 c where /3 op t is the optimal value of /? obtained through solving ((4T|) using such 
a linear approximation. In the second iteration, the resulting optimization problem has the same structure 
as the problem (14TI ) in which (3 C has to be set to obtained from the first iteration. This process 
continues and every iteration is obtained by replacing log(/3) at the iteration k by its linearization of type 
(l40l) around f3^ t found from the iteration k — 1. The POTDC algorithm for solving the problem ( f39b 
is summarized in Algorithm 1. 



Algorithm 1 The POTDC algorithm for solving the optimization problem (1391 ) 

Initialize: Select an arbitrary /3 C from theinterval [A m i n (Bj~ 1 B2), A max (B^ 1 B2)], set the counter k to 
be equal to 1 and choose an accuracy parameter e. 

while The difference between the values of the objective function in two consecutive iterations is larger 
than e. do 

Use the linearization of type d40l ) and solve the following optimization problem 

max log(tr(Ai • X)) + log(r) - t 
tr(Bi • X) = 1 
tr(A 2 • X) = r 
tr(B 2 • X) = p 

log(/3 c ) + l(p-p c )<t 

Pc 

X y 0. (42) 

to obtain xj } t , r^J, and /3 ( Q k J t . 
k = k + l 

Set X opt := X pj., and /3 C := f3^J t . 
end while 
Output: X opt . 



The following two lemmas regarding the proposed POTDC algorithm are of interest. First, the termi- 
nation condition in the POTDC algorithm is guaranteed to be satisfied due to the following lemma which 
states that by choosing /3 C in the above proposed manner, the optimal values of the objective function 
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of (|4TT) for X^pl-, r^pt, and are non-decreasing. 

Lemma 1: The optimal values of the objective function of the optimization problem ((4Tb obtained 
over the iterations of the POTDC algorithm are non-decreasing. 

Proof: Considering the linearized problem (ITil in the iteration k + 1, it is easy to verify that X^., 

(k) (k) 

Tq P {, and Pq P ( give a feasible point for this problem. Therefore, it can be concluded that the optimal 
value at the iteration k + 1 must be greater than or equal to the optimal value in the iteration k which 
completes the proof. ■ 
Second, it is guaranteed that the solution obtained using the POTDC algorithm is optimal due to the 
following lemma. 

Lemma 2: The solution obtained using the POTDC algorithm satisfies the Karush-Kuhn-Tucker (KKT) 
conditions. 

Proof: This lemma follows straightforwardly from a similar proposition in PTTl . ■ 
As soon as the solution of the relaxed problem 09l ) is found, the solution of the original problem (I35T ). 
which is equivalent to the solution of the sum-rate maximization problem ( f34T >. can be found using one 
of the existing methods for extracting a rank one solution. Among the existing methods are the ones 
based on solving the dual problem ll28l . which exploits the fact that the original problem (I35T ) with only 
two constraints is strictly feasible and has zero duality gap; the algebraic technique of ll30lk and the 
rank reduction-based technique of f29l which is also applicable for the problems with three constrains. 
Although the solution of (|39l ) is guaranteed to be optimal, it is still left to show that this solution is also 
globally optimal. 

V. An Upper-Bound for the Optimal Value 

Through extensive simulations we have observed that regardless of the initial value chosen for /3 C in 
the first iteration of the POTDC algorithm, the proposed iterative method always converges to the global 
optimum of the problem (|39l ). However, since the original problem is not convex, this can not be easily 
verified analytically. A comparison between the optimal value obtained by using the proposed iterative 
method and also the global optimal value can be, however, done by developing a tight upper-bound 
for the optimal value of the problem and comparing the solution to such an upper-bound. Thus, in this 
section, we find such an upper-bound for the optimal value of the optimization problem (|35l ). For this 
goal, we first consider the following lemma which gives an upper-bound for the optimal value of the 
variable (3 in the problem (I38T ). This lemma will further be used for obtaining the desired upper-bound 
for our problem. 
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Lemma 3: The optimal value of the variable j3 in (|39l ), denoted as /3 opt is upper-bounded by e^ q *~ p *\ 
where p* is the value of the objective function in the problem (l39l) corresponding to any arbitrary feasible 
point and q* is the solution of the following convex optimization problem^ 

q* = max log(tr(Ai • X)) + log(r) 

X,r,/3 

tr(Bi • X) = 1 
tr(A 2 • X) = r 

X y 0. (43) 



Proof: First note that since p* is the value of the objective function in the problem (1391 ) corresponding 
to an arbitrary feasible point, it must be less than or equal to the optimal value of problem d39l ). By fixing 
the variable f3 to /3 opt in the optimization problem (|39l ), the optimal value of the objective function does 
not change. Moreover, in the aforementioned case when (3 has been fixed to /3 pt> dropping the constraint 
tr(B2 • X) = /3 op t in that problem leads to the following optimization problem 

max log(tr(Ai • X)) + log(r) - log(/3 opt ) 

X,t,/3 

tr(Bi • X) = 1 
tr(A 2 • X) = r 

X y 0. (44) 



Noticing that the feasible set of the optimization problem (1391 ) is a subset of the feasible set of the newly 
introduced optimization problem (l44l . it is straightforward to conclude that the optimal value of the 
problem (1441 is bigger than or equal to the optimal value of the problem d39l and as a result it is greater 
than or equal to p*. Using (l43l . the optimal value of the optimization problem (l44l can be expressed as 
q* — log(/? op t) which is bigger than or equal to p* and, therefore, /3 opt < e^ q *~ p *^ which completes the 
proof. ■ 
Note that as mentioned earlier, p* is the objective value of the problem d39l that corresponds to an 
arbitrary feasible point. In order to obtain the tightest possible upper-bound for /3 op t, we choose p* to 
be the largest possible value that we already know. A suitable choice for p* is then the one which 
is obtained using the POTDC algorithm. In other words, we choose p* as the corresponding objective 
value of the problem (|39l at the optimal point which is resulted from the POTDC algorithm. Thus, 



4 Note that this optimization problem can be solved efficiently using numerical methods, for example, interior point methods. 
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we have obtained an upper-bound for /? opt which makes it further possible to develop an upper-bound 
for the optimal value of the optimization problem (l39l) . To this end, we consider the only non-convex 
constraint of this problem, i.e., log(/3) < t. Fig. [3] illustrates a subset of the feasible region corresponding 
to the non-convex constraint log(/3) < t where /3 m ; n equals A m i n (B ] ^ 1 B2), i.e., the smallest value of 
p for which the problem d39l ) is feasible, and /3 max is the upper-bound for the optimal value /3 opt 
given by Lemma 3. For obtaining an upper-bound for the optimal value of the problem ( f39b . we divide 
the interval [/3 m i n , /? max ] into N sections as it is shown in Fig. [3] Then, each section is considered 
separately. In each such section, the corresponding non-convex feasible set is replaced by its convex-hull 
and each corresponding optimization problem is solved separately as well. The maximum optimal value 
of such N convex optimization problems is then the upper-bound. Indeed, solving the resulting N convex 
optimization problems and choosing the maximum optimum value among them is equivalent to replacing 
the constraint log(/3) < t with the feasible set which is described by the region above the thin line in 
Fig. [3] The upper-bound becomes more and more accurate when the number of the intervals, i.e., N 
increases. 

VI. Semi-Algebraic Solution via Generalized Eigenvectors (RAGES) 

In this section we present RAGES as an alternative solution to the sum-rate maximization problem (l34l 
which is based on generalized eigenvectors. It requires a different parameterization than the one used in 
the POTDC algorithm and in some cases it is more efficient. 

A. Basic Approach: Generalized Eigenvectors 

To derive the link between (l34l and generalized eigenvectors we start with the necessary condition for 
optimality that the gradient of (l34l vanishes. Therefore, if we find all vectors g for which the gradient 
of the objective functions is zero, the global optimum must be one of them. By using the product rule 
and the chain rule of differentiation, the condition of zero gradient can be expressed as |[25l 

— — — ■ Ai ■ g + — — — • A2 • g 

Pn,i ■ Pn,2 Pn,i ■ Pn,2 

Pr,i ■ Pr,2 d . Pr,i-Pr,2 d 
= 62 5— • B l • S + 6 62- B 2 S- 

Rearranging (|45T ) we obtain 

(Ai + p sig • A 2 ) • g = S^i • (Bi + pnoi • B 2 ) ■ g (46) 
Pn,i 
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where p s i g and p no ; are denned via 

p R,i , Pn,i fAn . 

Psig = ~~. — and Pnoi = 75 — • (47) 
"fi,2 -nv,2 

It follows from (l46l ) that the optimal g must be a generalized eigenvector of the pair of matrices 
(Ai + p s ig • A2) and (Bi + p no \ • B2). Moreover, the corresponding generalized eigenvalue is given by 
Pr,i/Pn,i which is logarithmically proportional to the rate of the terminal one r\. Unfortunately, the 
matrices (Ai + p slg • A2) and (Bi + p noi • B2) contain the parameters p S i g and p noi which also depend 
on g and are hence not known in advance. Therefore, we still need to optimize over these two parameters. 
However, compared to the original problem of finding a complex-valued Mr x Mr matrix, optimizing 
over the two real-valued scalar parameters is already significantly simpler. The following subsections 
show how to simplify this 2-D search even further. 

B. Bounds on the parameters p s i g and p no i 

Since both parameters p s - lg and p no \ have a physical interpretation, the lower and upper-bounds for 
them can be easily found. Such bounds are useful since they limit the search space that has to be tested. 
For instance, p no \ can be expanded into 

Pn,i g H • Ji • g + Pn,i 
Pnoi = -~— = ~^r~i r~5 — ■ ( 48 ) 

PN,2 S • J 2 • g + Pn,2 

The quadratic forms can be bounded by using the fact that for any Hermitian matrix R we have 

A min (R) • ||g||| < g H ■ R • g < A max (R) • ||g||| (49) 

where A m i n (R) and A max (R) are the smallest and the largest eigenvalues of R, respectively. It follows 
from ((26ll that 

A m in(Ji) = and A max (Ji ) = A max (RAr jR ) • (af 5 ) 2 (50) 
where af^ is a short hand notation for ||h-^||2. Furthermore, in general the following inequity holds 

A max (R7V,i?) < Pn,r ■ M R (51) 

which for the case of white noise at the relay boils down to the following tighter condition A max (RAr j i?) = 
Pn,r- 

The relations (1501) and (IBTI ) can be used to bound (|48T ). Specifically, an upper-bound for p no \ can be 
found by upper-bounding the enumerator and lower-bounding the denominator, while the lower-bound 
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can be found by lower-bounding the enumerator and upper-bounding the denominator. This yields 

Pnoi < -5-*- ' Mr • (c4 ; ) • 7 Z + — — (52) 

"JV,2 "AT,2 

. (PN,R Ajr , (6)v2 2 , ^,2^ 1 /<Q s 

Pnoi > -5 M K -(a^) -7 +— — (53) 

\-nv,i Pn,iJ 

where = 1 1 S 1 1 2 and Mr can be dropped if the noise at the relay is white. However, an upper-bound 
for 7 2 is till needed. Due to the relay power constraint we have g^ ■ Q • g = Pt,r- Using the latter 
condition, the following bound can be derived 7 2 < Pt,r/ ^A m i n (Q). However, it is easy to check that 
this bound is very loose since for white noise at the relay we have A m i n (Q) = Pn,r and for arbitrary 
relay noise covariance matrices no lower-bound exists (the infimum over A m i n is zero). This bound is 
so loose because it is extremely pessimistic: it measures the norm of g in the case when only noise is 
amplified and no power is put on the eigenvalues related to the signals of interest. However, such a case 
is practically irrelevant since it corresponds to a sum-rate equal to zero. Therefore, we propose to replace 
7 2 in (EU) and (TJl bjS 

, 2 . = (»S /) ) 2 -(4 /) ) 2 -^ 1 -pt, 2 

(a[ f) y-P T>1 + (ai f) ) 2 -P Tt2 
In a similar manner, p s - lg can be bounded. In this case, the enumerator and the denominator have the 

additional terms Pt,2 ■ g H ■ K24 ■ g and Pt,i • g H ■ K-i )2 • g, respectively. A pessimistic (loose) bound is 
obtained by bounding these two terms independently, i.e., < g^ • K.2,1 • g < 7 2 • (a^) 2 ■ (af^) 2 and 
< g^ • K 1>2 • g < 7 2 • (a[ f) ) 2 • (4 b) ) 2 - This yields 

Psig < -5— • K ) 'K ) -7 +~B M R -{a\>) -7 + ^- L - (55) 

"JV,2 "iV,2 "JV,2 

. / Pt,\ , (/), 2 / (b)\2 2 , Pn,R , r , (6)%2 2 , ^V,2\ 1 

Psig> -5— -(ai ; ) 'K ) -7 +^ M R -(a^) - 7 +^ • (56) 

Again, these bounds are pessimistic since they assume that there exists an optimal relay strategy for 
which Pr i = P^,2 • (c4^) 2 ' ( a i^) 2 ' 7 2 but Pr i2 = 0, i.e., the rate of the second terminal is equal to 
zero. However, it is typically sum-rate optimal to have significantly more balanced rates between the two 
users. In fact, for the "symmetric" scenario when Pt : i = Pt,2, h$ = hj; > i = 1,2, and o^P = 
we always have Pr,\ = Pr >2 at the optimal point. Therefore, these bounds can be further tightened if a 
priori knowledge about the scenario is available. 

5 We have observed in all our simulations that this value poses indeed an upper-bound on the norm of the optimal solution 

gopt- 
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C. Efficient 2-D and 1-D Search 

Once the search space for p s - lg and p noi has been fixed, we can find the maximum via optimization over 
these two parameters using a 2-D search. In general, a 2-D exhaustive search can be computationally 
demanding, i.e., the complexity will be higher than that of the POTDC algorithm. However, as we 
show in the sequel, for the problem at hand, this search can be implemented efficiently. These efficient 
implementations are, however, heuristic since they rely on properties of the objective functions that are 
apparent by visual inspection. As we will see in simulations, the resulting RAGES algorithm performs 
as well as the rigorous POTDC algorithm in practice. 

Fig. [4] demonstrates a typical example of the sum-rate T\ + r 2 as a function of p s i g and p no ;. For 
this example we have chosen Mr = 6, Pr i = Pt,2 = Pt,r = 1, Pn,i = Pn,2 = Pn,r = 0.1 
and we have drawn the channel vectors from an uncorrelated Rayleigh fading distribution assuming 
reciprocity. By visual inspection, this sample objective function shows two interesting properties. First, 
it is a quasi-convex function with respect to the parameters p no \ and p s \ g which allows for efficient 
(quasi-convex) optimization tools for finding its maximum. Albeit this property is only demonstrated for 
one example here, it has been always present in our numerical evaluations even when largely varying 
all system parameters. Secondly, for every value of p S i g the corresponding maximization over p noi yields 
one maximal value which depends on p no \ only very weakly. This is illustrated by Fig. [5] which displays 
the relative change of the objective function v\ + r 2 for different choices of p QO \, each time optimizing 
it over p s \ g . The displayed values represent the relative decrease of the objective functions compared to 
the global optimum, i.e., for the worst choice of p no \, the achieved sum-rate is about 2 • 10 -5 = 0.002 % 
lower than for the best choice of p no \. Consequently, the 2-D search over p s - lg and p no \ can be replaced 
essentially without any loss by a 1-D search over p s i g only for one fixed value of p no \ (e.g., the geometric 
mean of the upper and the lower-bound). 

In addition, instead of performing the search directly over the original objective function n + r 2 , we 
can find an even simpler objective functions by using the physical meaning of our two search parameters. 
To this end, let us introduce a new parameter p s i g as a function of g as follows 

g H ■ Ai • g 

Mg) = g^-A 2 -g ' (57) 
Here g is the relay weight vector at the current search point (p s ig, Pnoi)- Then we know that in the optimal 
point g op t, we have /5 s i g (gopt) = Psig- This can be used to construct a new objective function 

^sig(Psig) Pnoi) — Psig(gopt) Psig (58) 
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Using the same data set as before, we display the corresponding shape of -Asig(Psig) Pnoi) in Fig- 
The red dashed line indicates the set of points where A s i g (p s i g , p no \) = 0. It can be observed that for 
every value of p noi , A sig (p sig , p noi ) is a monotonic function in p sig . Therefore, the bisection method can 
be used to find a zero crossing in p s j g which coincides with the sum-rate-optimal p s \ g for a given p ao \. 

D. Summary 

In summary, it can be concluded that the RAGES approach simplifies the optimization over a complex- 
valued Mr x M ft matrix into the optimization over two real- valued parameters which both have a physical 
interpretation. Even more, the 2-D search can be simplified into a 1-D search by fixing one of the 
parameters. The loss incurred to this step is typically small. In the example provided above, it is only 
0.002 %, but even varying the system parameters largely and using many random trials we never found 
a relative difference higher than a few percents. 

Moreover, the 1-D search can be efficiently implemented by exploiting the quasi-convexity of r\+r2 or 
the monotonicity of A s \ g (e.g., using the bisection method). Again, these properties are only demonstrated 
by examples but we have observed in all our simulations that the resulting algorithm yields a sum-rate very 
close to the optimum found by the exact solution and its upper-bound described before. This comparison 
is further illustrated in next section via numerical simulations. 

Comparing the POTDC and RAGES approaches, it is noticeable that the POTDC approach is absolutely 
rigorous, while the RAGES approach is at some points heuristic. The complexity of solving the proposed 
sum-rate maximization problem for two-way AF MIMO relaying using the POTDC algorithm is the 
same as the complexity of solving the semi-definite programming problem d39l and iterating over a 
single parameter /3. The typical number of iterations is 4-7. Alternatively, the complexity of solving 
the same problem using the RAGES approach is equivalent to the complexity of finding the dominant 
generalized eigenvector, which has to be performed for each combination of the parameters p s j g and p no \. 
Since, as has been shown, the search over one parameter only is sufficient, the complexity of the RAGES 
approach is typically lower than that of the POTDC algorithm, especially for the 1-D RAGES. 

VII. Simulation Results 

In this section, we evaluate the performance of the new proposed methods via numerical simulations. 
Consider a communication system consisting of two single-antenna terminals and an AF MIMO relay with 
Mr antenna elements. The communication between the terminals is bidirectional, i.e., it is performed 
based on the two-way relaying scheme. It is assumed that perfect channel knowledge is available at 
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the terminals and at the relay, while the terminals use only effective channels (scalars), but the relay 
needs full channel vectors. The relay estimates the corresponding channel coefficients between the relay 
antenna elements and the terminals based on the pilots which are transmitted from the terminals. Then 
based on these channel vectors, the relay computes the relay amplification matrix G and then uses it for 
forwarding the pilot signals to the terminals. After receiving the forwarded pilot signals from the relay 
via the effective channels, the terminals can estimate the effective channels using a suitable pilot-based 
channel estimation scheme, e.g., the LS. 

The noise powers of the relays and the terminals Pn,r, Pn,i and Pn,2 are assumed to be equal to 
a 2 . Uncorrected Rayleigh fading channels are considered and it is assumed that reciprocity holds, i.e., 
hjp = hf^ for i = 1, 2. The relay is assumed to be located on a line of unit length which connects the 
terminals to each other and the variances of the channel coefficients between terminal i, i = 1,2 and the 
relay antenna elements are all assumed to be proportional to 1/d", where € (0, 1) is the normalized 
distance between the relay and the terminal i and v is the path-loss exponent which is assumed to be 
equal to 3 throughout the simulations. [] For obtaining each simulated point, 100 independent simulation 
runs are used unless otherwise is specified. 

In order to design the relay amplification matrix G, five different methods are considered including the 
proposed POTDC, 2-D RAGES and 1-D RAGES algorithms, the algebraic norm-maximizing (ANOMAX) 
transmit strategy of ||32l and the discrete Fourier transform (DFT) method that chooses the relay precoding 
matrix as a scaled DFT matrix. Note that the ANOMAX strategy provides a closed-form solution for the 
problem. Also note that for the DFT method no channel knowledge is needed. Thus, the DFT method 
serves as a benchmark for evaluating the gain achieved by using channel knowledge. The upper-bound 
is also shown in all simulations. For obtaining the upper-bound, the interval [/3 m ; n , /3 max ] is divided in 
30 segments. In addition, the proposed techniques are compared to the SNR-balancing technique of ll33l 
for the relevant to the later technique scenario when multiple single-antenna relay nodes are used. 

A. Example 1: Symmetric Channel Conditions 

In our first example, we consider the case when the channels between the relay antenna elements and 
both terminals have the same channel quality. More specifically, it is assumed that the relay is located 
in the middle of the connecting line between the terminals and the transmit power of the terminals Pr \ 
and Pr t 2 and the total transmit power of the MIMO relay Pt,r are all assumed to be equal to 1. 

6 It is experimentally found that typically 2 < v < 6 (see 1341 p. 46-48] and references therein). However, v can be smaller 
than 2 when we have a wave-guide effect, i.e., indoors in corridors or in urban scenarios with narrow street canyons. 
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Fig. |7] shows the sum-rate achieved by different aforementioned methods versus a~ 2 for the case of 
Mr = 3. It can be seen in this figure that the performance of the proposed methods coincides with the 
upper-bound. Thus, the methods perform optimally in terms of providing the maximum sum-rate. The 
ANOMAX technique performs close to the optimal, while the DFT method gives a significantly lower 
sum-rate. 

B. Example 2: Asymmetric Channel Conditions 

In the second example, we consider the case when the channels between the relay antenna elements and 
the second terminal have better channel quality than the channels between the relay antenna elements and 
the first terminal and, and evaluate the effect of the relay location on the achievable sum-rate. Particularly, 
we consider the case when the distance between the relay and the second terminal, di, is less than or equal 
to the distance between the relay and the first terminal, d\. The total transmit power of the terminals, 
i.e., Pt,i and Pt,% and the total transmit power of the MIMO relay Pt,r all are assumed to be equal to 
1 and the noise power in the relays and the terminals all are assumed to be equal to 1. 

Fig. [8] shows the sum-rate achieved in this scenario by different aforementioned methods versus the 
distance between the relay and the second terminal denoted as 6%, for the case of Mr = 3. It can be 
seen in this figure that the proposed methods perform optimally, while the performance (sum-rate) of 
ANOMAX is slightly worse. 

As mentioned earlier, it is guaranteed that the POTDC algorithm converges to at least a local maximum 
of the sum-rate maximization problem. However, our extensive simulation results confirm that the POTDC 
algorithm converges to the global maximum of the problem in all simulation runs. Indeed, the performance 
of the POTDC algorithm coincides with the upper-bound. Moreover, the 2-D RAGES and 1-D RAGES 
are, in fact, globally optimal, too. The ANOMAX and DFT methods, however, do not achieve the 
maximum sum-rate. The loss in sum-rate related to the DFT method is quite significant while the loss in 
sum-rate related to the ANOMAX method grows from small in the case of symmetric channel conditions 
to significant in the case of asymmetric channel conditions. Although ANOMAX enjoys a closed-form 
solution and it is even applicable in the case when terminals have multiple antennas, it is not a good 
substitute for the proposed methods for the sum-rate maximization goal, because of this significant gap 
in performance in the asymmetric case. 
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C. Example 3: Effect of The Number of Relay Antenna Elements 

In this example, we consider the effect of the number of relay antenna elements Mr on the achievable 
sum-rate for the aforementioned methods. The powers assigned to the first and the second terminals as 
well as to the relay are all equal to 1. The relay is assumed to be located at the distance of 1/4 from 
the second user. Moreover, the noise powers at the terminals and at the relay antenna elements are all 
assumed to be equal to 1. For obtaining each simulated point in this simulation example, 200 independent 
simulation runs are used. 

Fig. [9] depicts the sum-rates achieved by different methods versus the number of relay antenna elements 
Mr. As it is expected, by increasing Mr (thus, increasing the number of degrees of freedom), the sum- 
rate increases. For the DFT method, the sum rate does not increase with an increase in the number of 
the relay antennas because of the lack of channel knowledge for this method. The proposed methods 
achieve higher sum-rate compared to ANOMAX. 

D. Example 4: Performance Comparison for the Scenario of Two-Way Relaying via Multiple Single- 
Antenna Relays 

In our last example, we compare the proposed methods with the SNR balancing-based approach of 
|[33l . The method of |33] is developed for a two-way relaying system which consists of two single- 
antenna terminals and multiple single-antenna relay nodes. Subject to the constraint on the total transmit 
power of the relay nodes and the terminals, the method of |[33l designs a beamforming vector for the 
relay nodes and the transmit powers of the terminals to maximize the minimum received SNR at the 
terminals. In order to make a fair comparison, we consider a diagonal structure for the relay amplification 
matrix G that corresponds to the special case of |[33l when multiple single-antenna nodes are used for 
relaying. It is worth mentioning that for imposing such a diagonal structure for the relay amplification 
matrix G in POTDC and RAGES, the vector gj\/2 xl = vec(G) is replaced with gA/ flX i = diag(G) 
and the matrices Aj and Bj,i = 1,2 are replaced with new square matrices Aj and Bj,i = 1,2 of 
size Mr x Mr such that A^ (m, n) = Aj ((m — 1) • Mr + m,(n — 1) • Mr + n) and Bj (m, n) = 
Bj ((m — 1) ■ Mr + m, (n — 1) • Mr + n), m,n = 1, ■ ■ ■ , Mr. Moreover, we assume fixed transmit 
powers at the terminals and choose them to be equal to 1. The total transmit power at the relay also 
equals 1 and the relay is assumed to lie in the middle of the terminals. Fig. [10] shows the corresponding 
performance of the different methods. From this figure it can be observed that the proposed methods 
demonstrate a significantly better performance compared to the method of |[33l as it may be expected. 
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VIII. Conclusions and Discussions 

We have shown that the sum-rate maximization problem in two-way AF MIMO relaying belongs to 
the class of DC programming problems. Although the typical approach for solving the DC programming 
problems is the branch-and-bound method, it does not have any polynomial time guarantees for its 
worst-case complexity. Therefore, we have developed in this paper two algorithms for finding the global 
maximum of the aforementioned problem with polynomial time worst-case complexity. The POTDC 
algorithm is based on a specific parameterization of the objective function, that is, the product of quadratic 
ratios, and then application of semi-definite programming (SDP) relaxation, linearization and iterative 
search over a single parameter. Its design is rigorous and is based on the recent advances in convex 
optimization. To the best of our knowledge, this is the first polynomial time algorithm for solving a class 
of DC programming problems rigorously. The RAGES algorithm is based on a different parameterization 
of the objective function and the generalized eigenvectors method, but may enjoy a lower computational 
complexity that makes it a valid alternative especially if 1-D search is used. The upper-bound for the 
solution of the problem is developed and it is demonstrated by simulations that both proposed method 
achieve the upper-bound and are, thus, globally optimal. 

The proposed POTDC algorithm represents a general optimization technique applicable for solving 
a wide class of DC programming problems. Essentially, the optimization problems consisting of the 
maximization/minimization of a product of quadratic ratios can be handled using the proposed POTDC 
approach. Moreover, the POTDC algorithm can be used for solving the optimization problems containing 
in any of the constraints a difference of two quadratic forms. Some relatively straightforward modifications 
may, however, be required. For example, if the problem is to optimize a product of more than two quadratic 
ratios under a single quadratic (power) constraint, the number of constraints in the corresponding DC 
programming problem will be more than three. Thus, the result used in this paper that even after relaxing 
the rank-one constraint in the step of SDP relaxation, it is possible to find algebraically an exact rank-one 
solution based on the solution of the relaxed problem, does not hold any longer. Then, randomization 
procedures will have to be adopted to recover a rank-one solution from the solution of the relaxed 
problem. In this case, such solutions obviously may not be exact, but all the results related to the SDP 
relaxation will apply. 

Other signal processing problems that can be addressed using the proposed POTDC approach are 
the general -rank robust adaptive beamformer with a positive semi-definite constraint [10;], the dynamic 
spectrum management for digital subscriber lines |[T2ll . the problems of finding the weighted sum-rate 
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point, the proportional-fairness operating point and the max-min optimal point for MISO interference 
channel lfT3l . the problem of robust beamforming design for AF relay networks with multiple relay 
nodes and so on. The extensions of the POTDC approach to some of the aforementioned problem is an 
issue of future research. 
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Fig. 1. Two-way relaying system model. 
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Fig. 2. Linear approximation of log(/3) around p c 
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Fig. 3. Feasible region of the constraint log(/3) < t and the convex hull in each sub-division. 
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Fig. 5. Relative change in sum-rate n + r-z versus p no [: optimizing over p s i g for every choice of p no i- The same data set as 
in Fig. [4] is used. 
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Fig. 6. Objective function A s i g (p s ig, Pnoi). The same data set as in Fig. 5] is used. The red dashed line indicates the points 
where A si g(p sig , p noi ) = 0. 
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Fig. 8. Sum-rate versus the distance between the relay and the second terminal d,2 for Mr = 3 antennas. The case of 
asymmetric channel conditions. 
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Fig. 10. Sum-rate versus a 2 for Mr = 3 antennas. The case of symmetric channel conditions. 
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